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Abstract 

Recent experiments on mesoscopic samples and theoretical considerations 
lead us to analyze multiply charged (n > 1) vortex solutions of the Ginzburg- 
Landau equations for arbitrary values of the Landau-Ginzburg parameter k. 
For n>l, they have a simple structure and a free energy T ~ n. In order to 
relate this behaviour to the classic Abrikosov result T ~ n 2 when k —* +oo, we 
consider the limit where both ii>1 and k> 1, and obtain a scaling function 
of the variable n/n that describes the cross-over between these two behaviours 
of T. It is then shown that a small-n expansion can also be performed and 
the first two terms of this expansion are calculated. Finally, large and small 
n expansions are given for recently computed phenomenological exponents 

characterizing the free energy growth with k of a giant vortex. 
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I. INTRODUCTION 



Recent experiments have demonstrated that vortices of charge n > 1, which are 

unstable in macroscopic type II superconductors ||, can exist in mesoscopic superconducting 
samples and can even be favored over configurations with multiple singly charged vortices. 
The Ginzburg-Landau description of superconductivity appears to be an adequate 

framework to analyze these results. The magnetization response of a mesoscopic sample can 
be analytically understood by adding a surface energy contribution to the Ginzburg-Landau 
free energy of a giant- vortex in an infinite system ||,[5). Moreover, full numerical solutions 
of the Ginzburg-Landau equations accurately reproduce the experimental findings flOHT2 



Motivated by these experimental and theoretical results, we analyze in the present paper 
'giant'(i.e. n > 1) vortex solutions of the Ginzburg-Landau equations. Well-known analyt- 
ical results have been obtained in the London limit when the Ginzburg-Landau parameter 
n — > oo |I3|,|7| and at the special dual-point value k = |7|JT^]. Here, we take advantage of 
the supplementary parameter n, the vortex charge, and provide a simple analysis of the giant 
vortices for arbitrary values of k. We begin in section 2, by recalling the Ginzburg-Landau 
equations and some elementary properties of their solutions. In section 3, we consider the 
large vorticity limit n ^> 1. The vortex structure takes the form of a circular normal core 
separated by a sharp boundary from the outside superconducting medium. As a conse- 
quence, the free energy T of a giant vortex is found to be proportional to its charge n in 
contrast to the classic Abrikosov's result T ~ n 2 in the London limit. In order to relate the 
two results, we consider in section 4, the double limit in which both the vorticity n and the 
Ginzburg-Landau parameter k are large and we obtain the scaling form T/2tc ~ nK$(«/n). 
The function $ provides an explicit interpolation between Abrikosov's result for k ^> n and 
the result of section 3 which is valid in the opposite limit n>K. To complete our analysis 
of giant vortex solutions, we consider in section 5, the somewhat more formal n — > limit, 
which has the advantage to be amenable to a simple and systematic expansion scheme; the 
free energy is obtained up to order n 2 . In section 6, we compare our perturbative results 
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(valid for arbitrary values of k) to the known exact results at the dual point k = l/y/2. 
Besides checking consistency, we obtain a perturbative expansion of the free energy in the 
vicinity of the dual point and the large-n expansion of phenomenological exponents intro- 
duced previously and computed numerically ||. The large n expansion proves to be fairly 
accurate for n values as small as 3 or 4. 



II. THE GINZBURG-LANDAU MODEL OF SUPERCONDUCTIVITY 

In the Ginzburg-Landau description of superconductivity the two unknown fields 
are the complex order parameter x/j = \ip\e tx and the potential vector A with V x A = B, 
where B is the local magnetic induction. These fields satisfy the following equations: 

2-7T 1 

-(v-z-a)V = t^(i-H 2 ), (i) 

Vx VxA=!|(^-4 (2) 

Here, the flux quantum 0o is given by 0o — hc/2e, and the two characteristic lengths A 
(penetration depth or London length) and £ (coherence length) appear as phenomenological 
parameters. The Ginzburg-Landau parameter k is defined as their ratio, k — j. We shall 
measure lengths in units of Av^2, the magnetic field in units of H c / \/2k = (po/AirX 2 and 
the vector potential in units of (f>o/2\f 7 2'K\. The Ginzburg-Landau free energy in units of 
if,?£ 2 /47r is then given by: 

Jr = ll B2 + K2{1 ~ ^ |2)2 + l( ^ " lA) ^ 2 - (3) 
A giant vortex of vorticity n is a solution of the Ginzburg-Landau equations with cylindrical 
symmetry ip(r,6) = f (r) exp(in6) , where (r, 9) are polar coordinates in the plane. The 
potential vector A can be chosen to lie in the plane and to have only a non-zero angular 
component Ag(r). It is convenient to introduce an auxiliary function g that represents the 
difference between the flux through a disk of radius r and the total flux, 

g(r) = rAg(r) — n i.e. B = (4) 

r dr 



The free energy can be rewritten in terms of / and g as 

^-jr{(D ,+ ^ + ^-^ + 5i(i)V- < 5 > 

The allied Ginzburg-Landau equations reduce to 

^ + -f-- 2 f9 2 = -2^f(l-f% (6) 

ar z r dr r z 

dr 2 r dr 

It is useful to note that a very simple form of the free energy H is obtained by an efficient 



S4$ = 2/V (7) 



use of and (|?p, 

—T = 2k 2 \ rdr(l-/ 2 ). (8) 

Z7T JO 

A vortex of charge n corresponds to functions / and g which satisfy /(0) = and 
g(0) = —n at the origin and which obey /(oo) = 1 and g(oo) = at infinity Linearization 
of Eqs. (H [?P around / = 1 and g = shows that there exists two exponentially growing 
and two exponentially decaying spatial modes at r = oo. In the same way, linearization 
for r close to zero shows that there is one diverging mode with / ~ r~ n and one neutral 
mode corresponding to changes in the vortex charge (note that at the level of Eqs. (§, 0), ra 
appears simply as a parameter and is not constrained to be an integer). Once one requires 
the diverging mode at r = to be absent and g(0) = —n, the expansion of / and g around 
r = depends on two arbitrary constants 

f(r)= (£) n + 0(r n+1 ), (9) 
g(r) = -n^-(0^+O(r 2n+2 ). (10) 

The length scales r and TZ are uniquely determined by the cancellation of the two diver- 
gent modes at r = oo they cannot be calculated from a local analysis near 0. Their 
determination requires the behaviors around r = and r = oo to be connected. This can 
be done numerically for arbitrary parameter values or analytically when n is either large or 
small as shown in the following sections. 



One can note from the definition of g (Eq. @) that ro is simply related to the value of 
the magnetic field at the position of the vortex: 



III. A GIANT VORTEX IN THE LARGE VORTICITY LIMIT 

We first consider the structure of a giant vortex of charge n for n ^> 1 and begin with 
simple estimates. It seems intuitively clear that the vortex core grows with its charge. So, for 
n ^> 1, one expects the vortex core to be much larger than the London penetration length. 
As a consequence, the magnetic induction should be approximately constant over the core. 
For a vortex of charge n and core size r c , one expects therefore \B\ = 2n/r 2 in the vortex 
core (the total flux divided by the core area with the chosen normalization) and B = 
outside. The magnetic energy of such a configuration is approximately / B 2 /2 ~ 2im 2 /r 2 . 
Correlatively, one expects ip = in the vortex core and ip = 1 outside and the corresponding 
energy contribution / k 2 (1 — IV'I 2 ) 2 — 2nK 2 r 2 . The total Ginzburg-Landau energy (j^) of such 
a configuration can therefore be estimated to be 

^ * -2 + + a ^ ^ 

where the last term on the right-hand side has been added to take into account the inter- 
facial energy of the transition layer between the vortex core and the superconducting bulk. 
Minimizing Eq. ([12]) with respect to r c gives the dominant order estimates for the vortex 
size r c and for its free energy T 

2 V2n 

r c ~ -— , (13) 

1 /2 

—F~nK\/2+ (^] <t(k) + 0(1). (14) 
2n \ k J 

The magnetic field inside the vortex core has the constant value 

B = V2k[1 + $L]. (15) 



That is, in the vortex core B is equal to H c plus a Gibbs-Thomson correction as expected 



at a curved normal /superconducting boundary (see e.g. |Iu , |r?|| ). In the next subsections, 
we first present numerical solutions of the Ginzburg-Landau equations for various values of 
n which confirm this simple picture of the giant vortex. This picture is then derived from a 
direct analysis of the Ginzburg-Landau equations in the large n limit. 



A. Numerical results 

Numerical solution of Eqs. (f||7]) consists in solving a two-point boundary value problem 
(at r = and r = oo) and can be achieved using several methods ||18|| . 

Shooting is generally easy to implement and robust when one mode needs to be can- 
celled at a boundary: a free parameter is adjusted (for example by dichotomy) at the other 
boundary until the required cancellation is obtained. In the present case, two diverging 
modes need to be cancelled at infinity which is more difficult to achieve. We kept a simple 
one parameter shooting by integrating from a large r = toward r = 0. Of the two free 
parameters at one was used to cancel the diverging mode at r = 0. The value of g at 
was then given as a function of the other parameter (which could be adjusted to obtain a 
particular value of n when desired). 

We also implemented a relaxation method. The system of ordinary differential equations 
is replaced by a set of finite difference equations satisfied on a mesh of points, and the bound- 
ary conditions just appear as equations satisfied by the points located at the extremities of 
the mesh. A multi-dimensional version of Newton's method provides a solution of these 
finite difference equations by an iterative procedure. This method requires the inversion of 
a matrix of size proportional to the number of points in the mesh. Because this matrix is 
block diagonal, it can be inverted in a very efficient manner. However, the efficiency of the 
relaxation method depends strongly on the starting point of the iteration: for example, it is 
useful to take the profile of a n-vortex as an initial guess for the (n + l)-vortex solution of 
the Ginzburg-Landau equations. 
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FIG. 1. a) The function / for k = 1.5 and n = 4,16,28,40,52. b) The different front profiles 
superposed on one another (/ = 1/2 is shifted to r = 0). 

Numerical solutions for different values of n are shown in Fig. 1. The plot of / shows a 
well-defined interface between a normal region (where / = 0) and a superconducting region 
(/ = 1) (see Fig. la). The interface width does not increase with n but remains of the order 
of 1/k. As expected, the interface has a well-defined limiting shape: in Fig. lb, the curves 
represented in Fig. la are shifted in such a way that they take the value 1/2 at the point 
r = 0. The shifted curves superpose on one another almost perfectly. Only the curve for 
the small value n = 4 shows a noticeable deviation from the limiting shape. 

The simple normal/superconducting interface picture for n ^> 1 is also confirmed by the 
graphs of g. In Fig. |^, B(r) = 1/r dg/dr is plotted for different values of n. The magnetic 
induction is close to H c = \/2k in the vortex interior and quickly decreases to zero in the 
interface region. Again, in the large n limit, the curves can be put over one another by the 
same shifts as those used to superpose the different /'s in Fig. lb. 
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FIG. 2. The magnetic field for k 2 = 1.5 and n = 4, 16, 28, 40, 52. 
B. Asymptotic analysis of the Ginzburg-Landau equations for n>l 



It is instructive to see how the large-n giant vortex structure arises in a direct analysis 
of the Ginzburg-Landau equations 



1. The vortex core 



Since g = —n at r = one expects g to be of order n in a whole region near r = 0. 
Therefore, at dominant order in n, Eq. (ffi) for / reduces to 



d 2 f 1 , 2 n 

dr r z 



(16) 



The solution of Eq. ( |l~6D can be obtained in a WKB manner but the important point is that 
/ is exponentially close to zero in the whole core region. Eq. (0) for g thus simplifies to 

d 2 g 1 dg 



dr 2 r dr 

At dominant order in n, g is therefore given by 



0. 



g(r) = -n [ 1 - ( — 



(17) 



(18) 



Although Eq. (|I8"D is identical to Eq. (|1~0|), it should be noted that it is not an expansion near 
zero but the full function g in the vortex core at dominant order in n. Of course, Eq. 
is nothing else but the constancy of B in the vortex core. 
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Outside the vortex core the medium is in the superconducting phase, / = 1 and g = 0. 
In order to match these different solutions of Eqs.(||[7|), we analyze the interface region (i.e. 
the boundary-layer in the usual matched asymptotics terminology). 

It should be noted that r c the vortex core size and tq which is obtained from the second 
derivative of g at the origin (Eq. (|?p) coincide at leading order in n. Their difference is of 
the order of the planar normal/superconducting interface width and depends on the precise 
criterion used to measure the vortex core size (i.e. on the precise definition of r c ) as shown 
below. 



2. Local equations near the interface 

The solution ([18]) vanishes for r = tq and therefore the assumption g ~ n under which it 
was derived breaks down at the interface. A different simplification of the Ginzburg-Landau 
equations can be made in the vicinity of r . We define a new variable x, centered around 
r , such that 

r = x + r . (19) 



We assume (and will verify at the end) that r ^> \x\. Substituting (|T9|) into (|],|7|) we obtain 



f l 2 +-^- 9 ^f = -^ 2 f(l-f 2 ), (20) 
dx A r ax rfi 

ax z tq ax 

Comparing the magnitude of the different terms in these simplified equations, we obtain 

/ ~ 1, g ~ ro, x ~ 1 and find that the first-order terms df /dx, dg/dx on the l.h.s. of ( PD||2TD 



are subdominant. Introducing the rescaled functions fo = f and Go = g/^o, one derives at 
dominant order in n (indicated by the subscript 0) the 'inner' equations which describe the 
behaviour of the order parameter and the vector-potential in the interface region 



d 2 fo 
dx 2 



2^/o(l-/o 2 ) + ^/ , (22) 



d 2 a 



of2 



2/oGo- (23) 



dx 2 
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As expected, these equations are identical to the Ginzburg-Landau equations in one dimen- 
sion (i.e. on an infinite line). Matching with the superconducting phase for x ^> 1 and 
the vortex core for x <C — 1 imposes the boundary conditions /o(+oo) = 1, Go(+oo) = 
and /o(— oo) = 0. This last condition implies, using (p3[) , that Go behaves linearly when 
x — > — oo 

G ~ ax + b. (24) 

In contrast to Eq. (H0), the local system (|22| , ^) is invariant by translation. Once this is 
fixed (for example by imposing the arbitrary criterion /o(0) = 1/2) the functions fo and Go 
and the constant b are uniquely determined (the value of a is independent of this arbitrary 
choice). Finding a and b appears to require an explicit solution. This can be avoided, at least 
for a, because the reduced system (^,^) retains the variational structure of the original 
equations and is invariant by translation in x. The locally conserved quantity £ associated 
with this continuous symmetry reads (that is, the magnetic field in the vortex core) |J 

Conservation of £ between — oo and +oo leads to: 

= £ (+oo) = E(-oo) = a 2 /A - k 2 /2. (26) 
Thus, a = \/2k and B = H c in the normal phase at dominant order in n, as expected. 

3. Matching the vortex core and the interface region 

In order to determine the constant and the interface position, we have to match the 
solutions in the vortex core and in the interface region. 
Rewriting fliTf ) in terms of the local variable x leads to 

I '-, /-, x \A 2n „ g In 

g ~ — n 1 — (1 H y ~ — x hence G = — ^ — (27) 

The family of solutions of (I22|j23|) is fo(x — x c ), Gq(x — x c ) with the asymptotic behavior 
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FIG. 3. a) The vortex-core size r c (dot-dashed curve) and the length tq (solid line curve) are 
plotted for k = 1.5 and for values of n ranging from 1 to 100. The dashed curve is the large-n 
asymptotic estimate (|l3|). The dots represent the next-to- leading order asymptotic expression 
(|38|), using the numerical value cr(1.5) ~ —1.172. b) Difference A ro between r$ and its asymptotic 
expression (|38|). 



G (x — x c ) = \/2k(x — x c ) + b 



(28) 



Identifying the two local expressions ( p7|) and ( p^) for the function g in their common region 
of validity r = r + x with \x\ -C r , and x — > — oo determines r and x c 
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The large-n asymptotic estimate of tq (which agrees with Eq. 
numerical results in Fig. p3 a. 



(29) 
(30) 

is compared with 



4- Large n limit of the energy 



The free energy of a giant vortex can now be calculated for n ^> 1. Using (|]) and 
integrating by parts, the free energy (||) becomes: 
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r<fr {T + K2 - K2 4 < 31 > 

Completing the square with the first two terms inside the brackets and using relation (|J) 
between B and g, leads to 



— jF = riKV / 2+ / rdr\ 
2-n Jo I 



The function 

is non-zero only in the vicinity of r . Hence, S can be considered to be a function of x = 
t — tq. Replacing the functions f(r) and g(r) in (^) by their dominant order approximations 



fo(x — x c ) and G (x — x c ), respectively, solutions of (|22j , |23|) , we obtain 

f°° f 00 f 00 | 1 / (IGq r-\ ^ I 

/ rS(r)dr — r S(x)dx = r / < — I — kv2 — k 2 /q > <ix = r cx(K). (34) 

io J ~oo J—oo I 2 V / 

The last integral is nothing else but the surface energy cr(«) of a one-dimensional domain 
wall between normal and superconducting phases |5, 19| , 20| (see appendix 0). Using (|34]) in 
(^2]) gives exactly the large n expression ( p!4[) of the free energy of a n giant- vortex with u{k) 
the surface energy of a one dimensional normal/superconducting interface. The expression 
(|i~4]) is valid for any value of k in the limit n — > oo. It is compared to the numerical results 
(for k — 1.5 ) in Fig. g. 



5. Finite-n correction to H c 

The leading fmite-n correction to -B in the vortex core can be derived from the previous 
results. Integrating (H) by parts leads to 

—f= K ' 2 [ 0O r 2 4)(r)dr with 0(r) = 2/(r)— . (35) 
27T jo ar 

The function <p(r) is localized at the position of the front and its width is of the order 1/k. 

Using the local variable x = r — ro defined in flT9|), the equation ( |55] ) can be rewritten, in 

the large n limit, as: 
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200.0 



FIG. 4. Free energy as a function of n for k = 1.5 (solid line). The horizontal dot-dashed 
line is the asymptotic value k\[2 ~ 2.121. The dashed curve is the large n behaviour (|i~4|), using 
the numerical value a (1.5) ~ —1.172. 



1 r°° 

-T = K 2 r\ + 2k 2 t q / x0(x)eb + O(l). (36) 

Z7T J-oo 

The integral in Eq. (p6|) can be expressed in terms of the interface energy (see appendix |A|): 
2k 2 j™J{x)xdx = 2k 2 jT ((1 - / 2 ) - = (37) 

Comparing ( |36"D with ( [14] ) provides the large n expansion of r up to the second order 



r 



n 



1/2 



K 



4/t 2 



1/2 / • 



(38) 



The vortex size r c , denned by the criterion f(r c ) = 1/2, can be meaningfully calculated to 
the same order: 



r c = r + x c 



(39) 



with x c given by Eq. (BO). 
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In order to verify numerically these relations, we first computed the value of the magnetic 
field at the origin using the following relation, that derives from (f|) and (^): 

,00 f2 g 



B(0) = -2 I" Lldr. 
Jo r 



This allowed us to obtain r , from (|ITD, with a good precision. The vortex size r c was found 
from the equation f(r c ) = 1/2. In Fig. 3a, the two numerical curves for r and r c are 
plotted as a function of n. To the leading order in n, these two curves differ by a constant 
x c ~ —0.797. . ., function of k only, as predicted by Eq. (|39|) . Solving independently the 
'inner' Ginzburg-Landau equations ( f22"U?3"D , we found that the unique value of b such that 
/o(0) = 1/2 is given by b ~ 1.691 . . . Hence the relation (|30"D between x c and b is verified. 
We also tested the accuracy of the large n asymptotic expansion of r to the first order (|i~3"D 
and to the second order (|38|). The difference between r and Eq. (p8|), hardly visible in 
Fig. 3a, is plotted in Fig. 3b. The agreement is very good. 



Substituting this value of r in (11), we obtain the field value inside the core 



B(0) = kV2 + ^L + o(z)- (40) 



nn \n 

Hence, when n — >• oo, the field at the origin is indeed the thermodynamic critical field 
H c = Ky/2 with a Gibbs- Thomson correction as expected (|15|) . 

IV. RELATION TO ABRIKOSOV'S FORMULA 

The self-energy of a unit vortex was first calculated by Abrikosov [0,0] in the large n 
limit 

^-T = log k + 0.081. (41) 

In the previous section we have studied the case n — > oo keeping k fixed and finite. In order 
to relate our results to the classical calculation of Abrikosov, we consider in this section the 
double limit n — >• oo and k —> oo keeping the ratio u = n/n finite and fixed. With these 
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FIG. 5. The scaling function $ as a function of the variable u (solid line). The dot-dashed curve 
and the dashed curve represent, respectively, the large u ( |65[ ) and small u (^) approximations. 



assumptions, we shall see that the Ginzburg-Landau equations decouple and that a scaling 
form is obtained for the free energy 



— T = nK$(«/n). 

Z7T 



(42) 



The function $ is plotted in Fig. [5]. 



A. The double limit n»l and k 1 



After rescaling the function g by a factor 1/n, and introducing the ratio u 
Ginzburg-Landau equations become: 

1 fd 2 f , ldf\ 2 1 2 a , 

- -2u f(l - ^^g - f 



n 2 \ dr 2 r dr / 
d 2 g 1 <ig 



2u 2 r 



dr 2 r dr ^ ^' 



ft/n, the 

(43) 
(44) 
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the boundary conditions are /(0) = and g(0) = — 1 at the origin, and /(oo) = 1 and 
g(oo) = at infinity. In the n —> oo limit, there are two domains where / is slowly varying 
and the derivative terms in Eq. (|4"5| ) can be neglected 

/ = for r <R, 

f=(l-?) ' for r>R. (45) 



As shown below, these two domains match through a boundary layer at r = R where / has 
a rapid variation but remains small. Matching the two slowly varying expressions (fi5"l) of / 
requires that f(R) = in the second one and gives 

g(R) = -V2uR. (46) 

Substituting ( fISD in fl44|) leads to the following closed equations for g: 

d 2 q 1 do . . 

for r < R —| ^ = 0, 47 

ar^ r ar 

for r > ^ _ 1^ = 2 (l- ^ \ (48) 
— dr 2 r dr \ 2u 2 r 2 j 



From (57), one has 



for r < R g(r) = -1 + (_j . (49) 

Continuity of the function (7 and of its derivative at i?, and equation ( ^6] ) imply that 

/ R\ 2 

g(R) = — 1+ ( — ) = -V2uR, (50) 
g'(R) = 2- 2 =\- 2V2u, (51) 



(the second equality of (|5"TD is derived using (5D)). The scaling limit for the free energy when 



n — > 00 and k — >■ 00 is obtained from (|) and (f45|) 
1 



2tt 



R 2 + 1 -[ 00 —dr) = nnu (r 2 + 2 T rh 2 (r) \ , (52) 

11 J Ft T J \ J R / 



where we have introduced a function h : 
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g = hV2ur. (53) 
This provides the explicit expression of the interpolating function $ flSq ) 

$(u) =u(r 2 + 2 I™ rh 2 (r)) . (54) 



>2 o / ..7.2, 

Ir 

In order to calculate explicitly $ one has to determine the position R of the front and the 
function h as a function of u — k/ti. In terms of h, Eq. ( [48|) becomes 

d 2 h 1 d/i 1 , , . , . 

— + -- -7^ = 2^(1-^), (55) 

ar z r ar r z 

and, from (^6|), the boundary conditions for /i are 

= -1 and h{oo) = 0. (56) 
Relation (^l|) implies, however, one more condition on h: 

The differential equation ( |55| ) with the two boundary conditions (^) and the supplementary 



condition (|57|) is overdetermined: there is a unique value of R such that all the conditions 
can be satisfied. Indeed, equation (|55| ) with the boundary conditions (|56| ) has a unique 
solution for any given value of R, i.e. for a given i?, the value of h'(R) is unique and can 
easily be computed numerically by solving (|55| ) by a shooting method; once h'(R) is known, 
u is determined as a function of i? from ( |5TD 

it 1 



(58) 



v/2 R + R 2 h'(R)' 

Inverting this relation gives R as a function of u (Fig. ||) and the function $ using (|54"1). 
This calculation can be performed analytically when u is either very large or very small as 
shown below. 

Before considering these limits, we complete the above analysis by giving the boundary- 
layer equation satisfied by / in the neighbourhood of r = R. It is convenient to introduce 
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FIG. 6. The inverse of R plotted as a function of the scaling variable u (solid line). The dashed 
curve represents the linear behaviour of 1 /R for large values of u (|6j) , and the dotted curve shows 
the asymptotic behaviour for small values of it (71). 

the local variable x = r — R. Assuming that / remains small in the neighbourhood of r = R 
and keeping the dominant contribution of each term of Eq. (pE5|), we obtain 

(59) 



L 2 {% + if] = 2u>f(-2h'(R)x + f). 
n z \ ax A Rax ) 



A consistent dominant balance between the different terms of Eq. (|59|) is obtained when 
f/(x 2 n 2 ) ~ fx ~ f 3 . Therefore, as assumed, / evolves on a short scale x ~ n~ 2 / 3 around 
r = R but remains small / ~ -nr 1 / 3 . With the rescaled coordinate and function, £ = 
x(4w 2 /i / (i?)) 1 / 3 n 2 / 3 , F = (n-u) 1 / 3 (v / 2/i / (-R)) _1 ^ 3 / ) tne boundary-layer equation reads 

d 2 F 



-£F + F 6 . 



(60) 



The two different functions of Eq. ( B5p can be matched through the solution of (the Painleve) 
Eq. P]). 
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B. The large u limit (k » n » 1) 



Recalling that R and r are positive, we see from flEJTp that R < \/{y/2u). Therefore, 
when u tends to infinity, R tends to 0. In the large u limit, equations fl4"T|) and ( f48|) reduce 
to: 



d 2 g 1 <ig 
<ir 2 r dr 



2g, 



(61) 



with the boundary conditions g(0) = —1 and g — > at infinity The solution of this equation 
is readily obtained: 



9{r) = -r 



\/2ifi(r\/2) and therefore /i(r) = -K x {r\f2)ju. 



(62) 



where i^i is the modified Bessel function of order 1. Using the behaviour of K\ in the 
vicinity of zero (see Appendix B), we deduce that 

1 



h'(R) 



V2uR 2 ' 



(63) 



Substituting 



in fl58T) leads to 



R 



\/2< 



u 



(64) 



Inserting (^2j) and (Q) in the the expression ( |o*^ ) for the free energy, we obtain 



1 ^ riK 
—T ~ — 

2tt u 



l/u 



r(Ki{r)) 2 dr 



Tin I 1 
— 1 + 



2u 



/v r, ( - ) I\ > (- 



,u 



~ — (lnu + ln2 - 7) ~ n 2 (In I- ) + 0.116 ) , 



(65) 



where the K u are modified Bessel functions and 7 is the Euler constant |[21|| . This relation 
generalizes the classical result (|4l| ) of Abrikosov to the case of a giant vortex. Indeed, the 
structure of the giant vortex is very similar to the n = 1 case with / varying on a fast scale 
of order ti/k and g on a slow scale of order one. However, the structure of / is simpler in the 
large n limit (being for r less than nj[\plH) and linked to g otherwise) than for general n. 
This is the reason why no new constant needs to be numerically determined in (|65|). 
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C. The small u limit (n > k > 1) 



In the small u case, the length R tends to infinity. Defining a local variable x = r — R, 

the differential equation fl55l) for h becomes 

d 2 h 1 dh 1 
dx 2 R + x dx (R + x) 2 

with the following boundary conditions 



+ iTIZ " 7777"^ = 2M1 - ^), (66) 



h(x = 0) = -1 and fc(oo) = 0. (67) 



When R — > oo, the equation (|66|) reduces to: 



2/i(l - /i 2 ). (68) 



dx 2 

The solution that satisfies the boundary conditions is: 

1 /2 I 

h(x) = -V2 (l - tanh 2 (v / 2(x + x ))) with tanh(v / 2x ) = (69) 
^From this expression, we derive that 

h'{r = R) = h'(x = 0) = 1. (70) 
Substituting this result in (|58| ) leads to the small u behaviour of R 

The subleading behaviour of R as a function of u is found by retaining only the first order 
term in 1/R in (|66|): 

d 2 h 1 dh . , 9 . . „ s 

+ = 2,1(1 < 72) 



with the same boundary conditions as in (|6y|). This equation can be solved perturbatively, 
but we only need to determine the correction to equation ([701); i.e. we need to calculate the 
coefficient C such that h'(r — R) = 1 + C/R. We use an argument of energy conservation. 
Defining 
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we readily obtain from ([F2|) that 

dS 2 1 (dh\ 2 
dx R \dx J 

Integrating this relation from to oo and evaluating the energy £2 at both ends leads to 

00 / dh x 2 



°=k U dx - (73) 



Using the Oth-order solution ( j69|) for h in ( |73|) we obtain 



h'(r = R) = h!{x = 0) = 1 + ( 2v ^ o - ) 4 (74) 



R 



Substituting this result in (p§) and solving for i? leads to 



R ^ V2 (l - 2 -V4^tl u l/2 + . (75) 



Substituting the expressions fl75|) and (^) for i? and /i, respectively, in (|52|) we obtain the 
free energy 

^-T - nn (uR 2 + 2R h 2 (x)dx ) j ~^ v / 2^1-^(v / 2-1)2 1/4 ^) 1/2 + C(m)^ . (76) 

These results are consistent with those obtained in section 3. In fact the small u limit 
amounts to first taking the limit n —>■ 00 and then taking k — > 00. For instance, it can be 
shown that equation (|38|) reduces to ([75]) in the k — > 00 limit. Moreover, comparing (ff6|) 
with (|14|) we retrieve the asymptotic behaviour of the surface tension in the large k limit 

MM- 



-(v 7 ^ - l)/t 2 ~ -0.7810k 2 . (77) 



V. A VORTEX IN THE SMALL VORTICITY LIMIT 

The winding number n must be a priori an integer because the phase of the wave function 
ip has to be a 27r periodic function. However, as already noted, in the cylindrically symmetric 
Ginzburg-Landau equations (P,0), n appears only as a parameter in the boundary condition 
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at the origin and can be given any real value. In this section, we take n to be close to and 
determine perturbatively the solutions of (^J^) and the free energy. 

Away from the origin r = 0, / is close to one and g is small. Thus, we expand / and g 
as follows, 

/ = 1 + nfi + n 2 f 2 + ... and g = ng 1 + n 2 g 2 + ... (78) 

The functions fi,gx, f 2 , g?,--- satisfy a hierarchical system of linear differential equations that 
can be solved recursively by imposing the boundary conditions: fi—>0 and gi — ► when 
r — > oo, for alH > 1. 

Thus the first order terms j\ and g\ are solutions of 

d 2 fi , ld/ x 



+ z^r = 4 « 2 A. ( 79 ) 

2g x . (80) 



(ir 2 r dr 
d 2 g 1 1 



<ir 2 r dr 

Implementing the boundary and the matching conditions, we obtain 

h = AK {2kt) , 9l = BV^rK^r), (81) 
where A and B are two undetermined constants at this stage. 



The expansion (|78f) breaks down in the vicinity of the origin since /(0) = 0. The 



expansion of / near r = is given by (§) 

/W= (^)" + 0(r" +1 ). (82) 
The two expressions (BT]) and fl8"2| ) can be matched in the parameter range where r < 1 



but n| log(r)| 1. When r < 1, the correcting terms in ( 82Q are negligible and when in 
addition n| log(r)| 1, one obtains 

/ = 1 + nln^ + 0(n 2 ,r). (83) 



In the same parameter range, expansion of Eq. ( jST|) gives 

/ = 1 + nA (-ln(«r) -7 + C(r 2 lnr)) . (84) 
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Consistency of these two expressions requires A = — 1 and 

^ 1 , x 0.561 

K = - exp(-7) ~ , (85) 

K K 

7 being the Euler constant. 

We can proceed in a similar way for the function g. Near r = 0, the expansion of g is 

- 1 + y)-^rryy + ° (rWy (86) 

For r < 1 and n\ log(r)| 1, this gives 

(1 77, 77 7^ \ 

~2 + ^|+ 2 -^ln^ + C(n 2 ,r)J . (87) 

Again, an alternative expression is obtained from the small r behavior of (|8~ID 

g = nB ^l + r 2 lnr + y(27-l-ln2) + 0(r 2 lnr)^J . (88) 

Comparing Eqs. (|8"?D and ( |88| ) gives B = —1 and 

r 2 = 2n -2n 2 ln(2/t 2 ). (89) 

Having determined the small-n expression for /, the free energy limit can be calculated 
from (H). Actually, the expression ( jSl]) for / is sufficient for this purpose (the range where it 
is not valid gives an exponentially small contribution) and by performing the integral over 
modified Bessel functions, one obtains J r /(27r) = n + 0(n 2 ). 

A similar procedure can be carried out for the second order terms. The expressions for 
/2 and gi are more intricate and are given in Appendix B. Using Eq. (|8D, they provide the 
free energy in the small n limit up to order n 2 

-Ljr = n + n z \ g( K V2) + 0(n 3 ) . (90) 

27T 
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VI. THE DUAL POINT OF THE GINZBURG-LANDAU EQUATIONS 



A. The dual point equations 

Eq. (pOD suggests that the value k = l/v2 has special properties (for instance, for 
k = l/y/2 the functions fi and g\ of ( [81] ) satisfy simple identities such as dfi/dr = —g±/r). 
It is indeed well-known that at this 'dual' point the second order Ginzburg-Landau equations 
reduce to first order equations, leading to special relations between the functions fi and 
Moreover, the free energy at the dual point has been calculated exactly and is identical to 
the topological number n. We recall briefly these special properties. 

At the dual point the free energy ([5]) can be written as follows: 

In Jo | (dr r ) 2 ^ ^ r dr) | Jo {c/r^ ^ ^} 

=f{(!^) 2 ^( (1 - /2) -^)> dr + "' (91) 

The minimal free energy is obtained when the following first order system is satisfied (Bo- 
gomol'nyi equations |0,|22],[14|] ) : 



Substituting (03) and (|53J) in ( [51] ) we obtain the minimal energy 



2vr 

Equation fl93|) implies that 



— ^ = 7i. (94) 



g = —n H / rf 2 (r)dr i.e. = 2n. (95) 

2 Jo 

Another remarkable consequence of the Bogomol'nyi equations is that the surface energy 
<t(k), defined in (p?4Tj, vanishes identically and changes its sign at the dual point (this is the 
reason why the dual point separates Type I from Type II superconductors 0): 
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cf{k 



1 



) = 0. 



(96) 



V2 



All the calculations that we have carried out in the previous sections are consistent with 
these (non-perturbative) properties of the dual point. In the small n case of section 5, 
one can verify, using (jSll) and the expressions given in the Appendix, that the equations 
( |92l , |93|) are satisfied order by order by the expansions ( |78|) of / and g at the dual point. 
Moreover, the n 2 correction to the free energy in ( p0[) and to r in (|89f) vanishes identically 
at k = 1/ y/2, as expected. In the large n case, the expression ([14]) for the free energy reduces 
to fl9~4l) because the n 1//2 correction disappears thanks to the vanishing of the surface energy; 
for the same reason, the formula (|3^) simplifies to r = (2k) 1 / 2 . 

At the dual point, the vortices do not interact and the free energy does not depend 
upon their location. This makes it possible to obtain exact n-vortices solutions, in the large 
n limit, for arbitrary locations of the vortex cores [pi. However, when k deviates from 
the special value l/y/2, vortices start interacting and their interaction energy is extremal 
when all the vortices are located at the same point. We now evaluate the energy of such a 
configuration with cylindrical symmetry when k is close to the dual point. 



experiments on mesoscopic superconductors, values of k are close to the dual point value 



useful to study analytically the magnetization of a sample as a function of the applied field. 
In the vicinity of the dual point, a relevant quantity is the relative growth of the free energy: 



B. Free energy of a giant vortex in the vicinity of the dual point 



We apply our results to the case where k is close to the dual point value k 



l/y/2. In 



|]J9|] and in previous theoretical studies the special case k 



1 / y/2 was found to be very 




when k 




1 



(97) 



This relation is a local form of an empirical scaling of the free energy 




(98) 
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found in ||, where the exponents a(n) were computed numerically. Differentiating the 
Ginzburg-Landau free energy ([5]) with respect to the parameter k at the dual point we 
obtain: 

1 t2\2 j I [°° dr { dg x 2 



a(n) = -J o d-/Vr* = -y o - ^ . (99) 



The last relation was derived from the Bogomol'nyi equation (|93|) . 

We now explain how a large n expansion of the exponents a(n) can be derived from our 
previous results. Integrating Eq. ( p9[) by parts, we find 

1 roo r % (if 

a(n) = - —ip(r)dr with ip = 4/(1 - f 2 )^-. (100) 

Tli J (XT* 

The function ip(r) is localized at the position of the front and its width is of the order 1/k. 
Hence we can expand (|100|) using the local variable x = r — and in the large n limit 
replace / by f (as in section 3.5): 

= ~ T 4 + - f°° ^)dx with V = 4/ (l - / 2 )^. (101) 
Integrating the last term by parts, we obtain: 

/0 /*oo 
- /o 2 ) 2 - l} + J dx(l-f$) 2 

= /_"^{(l-^) 2 -^} = -/_"^(l-^- ( 102 ) 



oo 

.ri 



In the last equality we used the fact that the functions /o and Go satisfy the following 
Bogomol'nyi equations at the dual point |TB|] : 



^ = (1-/ 2 ) and ^ + G /o = 0. (103) 

The function t> defined by / = exp(— v/2), satisfies a second order differential equation that 
can be solved explicitely []16[: 



- — r = 1 — exp(— 1>) which implies — = — 2 fe u — 1 +v) ^ . (104) 

2 dor ctx v 7 

Now changing the variable x to we rewrite the last integral in (|102|) as 
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_ roo e -* (1 e -v )d = ^ .oc _ : y/ 2 ^ _ 548 2. (105) 

Substituting this result in ( |101| ) and knowing that r 2 = 2n exactly fl9~5|), we find an asymp- 
totic expansion for a(n) in the large n limit: 

a(n) = i-^%B + o(±). (106) 



n 1 / 2 \n, 

We have not calculated the subleading behaviour of the free energy but our numerical results 
provide an estimate of the higher order term in ( |106| ): 

a(n) = l-^^ + - + o(-^) with K~ 0.125. (107) 

From (|106|) and (|97|), the expansion of the free energy near k = is found. This allows us 
to retrieve, using (|14D, the local behaviour of the surface energy in the vicinity of the dual 
point, first calculated by Dorsey fL6|: 



a(K) ~ —IK 



V2-l)J o dv (e~ v - 1 + vj e v ~ -0.5482(/t V / 2 - 1). (108) 



In the same way, from the small n expansion of the free energy up to the second order, 
we derive that (see Appendix B): 

a(n) = n-4n 2 r 2 ^(r)Ki(r)o?r + 0{n 3 ) ~ n - 1.5626n 2 + 0{n 3 ). (109) 
Jo 

We have computed numerically from (^7|) the function a(n) for n ranging from 1 to 250 
and some significant values are given in Table 1. The numerically computed and the large 
n asymptotic expansion ( |106|) of a(n) agree to better than 7% even for n as small as 4. 
In particular we notice from ( |106|) that a(n) — > 1 when n — > oo. In the n — > limit, the 
second order expansion ( |109|) provides a fairly good approximation for n < 0.25. In Fig. 7, 
numerically computed values of a are compared to large n and small n expansions. 



VII. CONCLUSION 



We have analyzed the structure of a giant vortex of winding number n in an infinite plane 
for arbitrary values of the parameter k in contrast to previous analytical results obtained 
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n 


a(n) 


n 


a(n) 


n 


a(n) 


1 


0.415 


12 


0.789 


50 


0.893 


2 


0.542 


13 


0.797 


60 


0.902 


3 


0.611 


14 


0.804 


80 


0.915 


4 


0.655 


15 


0.810 


100 


0.924 


5 


0.687 


16 


0.816 


120 


0.930 


6 


0.711 


17 


0.821 


140 


0.935 


7 


0.730 


18 


0.826 


160 


0.940 


8 


0.746 


19 


0.830 


180 


0.943 


9 


0.759 


20 


0.834 


200 


0.946 


10 


0.770 


30 


0.863 


225 


0.949 


11 


0.780 


40 


0.881 


250 


0.951 



TABLE I. Numerical values of the function a{n) obtained from (p9|) for different n in the range 



1 to 250 
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S 0.6 




200.0 



FIG. 7. The function a(n) is plotted for n between 1 and 200, and in inset for values of n 
between and 1 (solid line curves). The dashed lines represent the asymptotic expansions for large 



n (|i06D and small n (|i09|) 



only in the London limit or at the dual point. The vortex and magnetic field profiles are 
computed by solving the Ginzburg-Landau equations which minimize the free energy of the 
system. These numerical solutions are compared with analytical results derived in the cases 
where the vortex multiplicity n is either very large or very small. In particular, a simple 
structure has been found for large n, its relation to the classic result of Abrikosov has been 
elucidated and the perturbative expansions are found to agree well with previous numerical 
computations of a giant vortex free energy. We hope that some of these results will prove 
useful for the current very active experimental investigations of mesoscopic superconductors. 

Acknowlegments. This work was initiated when two of us (V.H. and K.M.) were 
attending the workshop "Ginzburg-Landau models, vortices and low-temperature physics" 
at the Lorentz Institut, Leiden. It is a pleasure to thank the organizers, of this conference, 
particularly Eric Akkermans, for stimulating discussions and the Lorentz Institut for its 
hospitality. K.M. would like to thank S. Mallick for his constant help during the preparation 



29 



of this work. 



APPENDIX A: EXPRESSIONS OF THE NORMAL/SUPERCONDUCTING 

INTERFACE ENERGY 

The energy per unit area of a planar normal/superconducting interface is obtained from 
the Ginzburg-Landau energy (H) under the form [|] 

^ - c { (f ) 2 + ^ + ^ - ^ + H§) 2 - } < ai » 

where fo and Go satisfy Eqs. (|^. |23|). We notice that for a magnetic field oriented along 



the z-direction, Go represents the ^/-component of the potential vector. The expression (|A1|) 
can be written in different forms by making use of the identities: 

II (f = II { 2 " 2/ ° 2(1 " /q2) ~ Glf ^ dX ' (A2) 
dx ([ ^ 2 - m/2^| | = -2 jT Gg/'ds, (A3) 



The first two identities are obtained (e.g. [|2D[) from Eqs. (^,^) by integration by parts. 
The third identity results from the conservation law (E3). Using (|A2|) to eliminate / (dfo/dx) 2 



from Eq. ( |A1| ) provides the expression ( |5^ ) of the main text. Similarly, using Eqs. 
to eliminate f (dfo/dx) 2 , J(dGo/dx) 2 and / /q^o f rom Eq. (Al) leads to the alternative 
expression (|37D of the main text: 

4 .2 f°° (m ^ 1 d Go^ 



^»=3 K /-4 (1 - /o) -^r' (A5) 



Finally, we note that the integral of x <p in (j37| ) is transformed into Eq. ( |A5| ) by integrating 
separately by parts the integrals between — oo and and between and +oo: 

4>(x)xdx = [ ( (x + Xc)-jr- ) dx = ( (l — /g) dx — / f^dx 

M- fS) -^0) dx (A6) 
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(the last equality is obtained by adding under the complicated form 

= / dx[ — — - — j + / dx 1 p-T^ • A7 

APPENDIX B: SECOND ORDER COMPUTATION IN THE N— > LIMIT AND 

SOME USEFUL FORMULAS 

We first recall some useful asymptotic formulae for the modified Bessel functions: 

K (r) = - In 7 - - 7 + 0(r 2 \nr), 

K x {r) = l - + T - In T - + ^(2 7 - 1) + 0(r 2 lnr), 

K 2 (r) = ^ + 0(r- 1 ). (Bl) 

— — = —K\(r) and — (rXi(r)) = —rK (r). 
dr dr 

^ = h{r) and -^-(r/^r)) = r/ (r). (B2) 
or ar 

We now explain how the matching procedure is carried out up to the second order in the 
case n — > 0. The functions fa and 52 solve the linear system: 

d 2 fa 1 dfa 2 # 2 2 2 
-tt + -^- = 4k / 2 + -a +6k /i, 
ar^ r ar r z 

^2 1^2 

ar^ r ar 

Solving this system by 'variation of constants', we obtain: 

fa = AKq{2kt) + BI q (2kt) - J (2«r) J™ uduK (2nu) (^K^uV^) 2 + 6k 2 K (2ku) 2 ) 

- K (2nr) J* uduI (2Ku) (2K 1 (m v / 2) 2 + 6/t 2 fsr (2/™) 2 ) , (B3) 

^ 2 = CK (2Kr) + J D/ (2/tr)-4v / 2r/ 1 (v / 2r) /°° tf (2Ku)(iiri(>/2u)) 2 udu 

-4v / 2rX 1 (v / 2r) f K Q (2Ku)K 1 (V2u)h(\f2u)udu. (B4) 

JO 
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Since f'2 and #2 must be finite at infinity, we have B = D = 0. The two other coefficients 
are found by matching with the inner expansions in the vicinity of zero. For this purpose, 
we need the local expansions of / and g in the vicinity of up to the order n 2 . Using ( |82| ) 
and (HH), we obtain 

f{r) = l + n(\n^ + 0(rlnr)\ + y ((In ^) 2 + 0(r lnr)) + C(n 3 ). (B5) 



g (r) = -n[l + 0(r 2 \nr)) +n 2 (r 2 ln— (1 - In—) + 0(r 2 ) ) + 0{n 6 ). (B6) 



Because the coefficient of n 2 in the local expansion ( P36|) of g tends to zero when r vanishes, 
one must have (72 — ^ when r — > and therefore C = 0. To obtain the coefficient A, we must 
determine the divergences of f'2 given by (|B3|) when r — > 0: the term AKq{2kt) produces 
a singular term — Alnr when r — > 0; from the expansions of the Bessel functions ( |B1|) , we 
find the diverging part of the first integral in ( |B3|) to be — J r °° 2uduK (2Ku) Ki{u\/2) 2 which 
is equivalent to j 2nr (In | + rfdu/u; thus the first integral gives a singular term equal to 

( lnr ) 2 / M 
h (—7 — uxk) Inr ; 

similarly the singular term due to the second integral in (|B3|) is given by: 



(lnr) 2 - lnr (-7 - hi ft + 2 f 1 wdn (j (2K«)lifi(W2) 2 L 

Hence the total diverging part of f2 is: 

- lnr (A + 2 jT 1 udu (/ (2km)ATi (wv^) 2 ~ ^2)) ■ (B7) 

Identifying this expression with the singular part of the n 2 term in the expansion ( |B5[ ) of / 
and using fl55|) provides the value of A and leads to the second order terms in the small n 
expansion of / and g: 

/•oo 

Iku) 2 



uduK (2Ku) (2#i(W2) 2 + 6K 2 K (2t 
-K (2Kr) J\duI (2Ku) (^K^uV^) 2 + 6K 2 K {2t 

with A = -7-I11/C- 2 / «d« (j (2«;u)tfi(W2) 2 - At) ; (B8) 

jo V 2u 2 J 



IKU) 2 
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roc 

g 2 = -4 v / 2r/ 1 (v / 2r) / K {2nu){K 1 {V2u)) 2 udu 

Jr 

-iV^rKUV^r) f K (2Ku)K 1 (V2u)I 1 (V2u)udu. (B9) 
Jo 

We can now derive the small n asymptotic behaviour of a(n). From the explicit relation 
]), we obtain, up to the second order in n: 

1 f°° dr ( dg\ 2 f°° dr ( dgi\ 2 2 r°° dr f dgi dg 2 \ , 3> 



^ ^ n Jo r \dr J Jo r ( dr J Jo r \dr dr J ^ ^ ^ ^ 

where the functions gi are calculated at the dual point k — 1/ \/2. Recalling that g\ = 
-V2rK x (^/2r) (see |81|) and using ( |B2|) we have 



dr ( dgi 



2 



oo 



n Vl^l =2 /" xK Z( X ) dx = 1 - ( B11 ) 







From (B2) and (]B9| ) , we obtain dg±/dr = 2rKo(^/2r) and 

^ = 4r\-I (V2r) f°° K (u)(K 1 (u)) 2 udu + K (V2r) [ * Koi^K^h^udu) . 
dr [ Jy^r Jo J 

(B12) 

The coefficient of n 2 in ( B10|) is therefore given by: 

2 r - (^^] = A T drK (V2r)^ = -8 f° dvvK (v)I (v) l°° K^K^fudu 
Jo r \dr dr J Jo dr Jo Jv 

roc rv 

+ 8 / dvvK 2 {v) / K {u)K x (u)Ii{u)udu. (B13) 

JO JO 

From ( |B2| ) we can verify that 

2 2 

j dvvKl{v) = V -{Kl-Kl) and | dvvK {v)I {v) = j{K I + K x h). (B14) 
Integrating ( |B13| ) by parts leads to: 

/"OO /*00 /"OO 

-8/ dvvK (v)I (v) K (u)(Ki(u)) 2 udu — —4 / dvv 3 (K I + K^KqK 2 , 

Jo Jv Jo 

/•OO /"f /"OO 

8/ dvvKhv) K (u)K 1 (u)I 1 (u)udu = -4 dvv 3 (K 2 - KhK^h. (B15) 

JO JO JO 

Summing these two expressions and using Koh + K\Iq = 1/v leads to 

2 l X ± { ^f) = -4 f ivv*Kl {v )KM- (B16) 
Jo r y dr dr J Jo 

This concludes the derivation of (|109|). The free energy in the small n limit up to order n 2 



is derived by similar calculations. 
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